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FOREWORD 


This final technical report was prepared by Northrop Norair, Hawthorne, Cali- 
fornia on Contract AF 33(657)-8851 for the Aeronautical Research Laboratories, 
Office of Aerospace Research, United States Air Force, The work reported herein was 
accomplished on Task 7065-01, "Fluid Dynamics Facilities Research" of Project 7065, 
"Aerospace Simulation Techniques Research" under the technical cognizance of 
Lt Arthur Wennerstrom of the Fluid Dynamics Facilities Laboratory of ARL. Report 
number NOR 63-59 has been assigned to this report for internal control, 
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ABSTRACT 


This report describes an IBM FORTRAN program for calculating supersonic flows 
through curved channels. Included in the program are: 1) a second order 
Belotserkovskii procedure for calculating the flow field around the forward part 
of blunt blades, 2) a general method of characteristics procedure for calculating 
supersonic flow fields, 3) automatic procedures for calculating interactions of 
shocks of opposite families, shock interactions with a wall and with a vortex 
sheet or slipstream, 4) a finite difference laminar boundary-layer procedure for 
calculating viscous region of the flow. Three sample calculations as well as the 
listing of the FORTRAN program are presented. 
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(BLUNT BODY SOLUTIONS) 


SYMBOL DESCRIPTION 
x, y Cartesian coordinates 
r, 6 Polar coordinates 
M Mach number 
Ww, P,P Dimensionless speed, pressure, density 
Wy Wy Velocity components (x, y directions) 
u, V Velocity components (r, 6 directions) 
e=p/p * Entropy function 
y Stream function 
y Adiabatic index 
k (y-1)/2y 
0 Angle between tangent to shock wave and incident stream 
€(6) Distance between body and shock measured along radius vector 
t TV 
1 
r aw) 
8 puv 
g 2kp + pv 
H kp + pue 
h Tu 


Other symbols are defined in the text. Velocities are made dimensionless in 


terms of the maximum velocity, pressure and density in terms of corresponding 
stagnation values. 
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SYMBOL 
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Pp - Local Static Pressure 
R - Local Total Pressure 
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I_ INTRODUCTION 


The behavior of supersonic flow through cascade depends strongly on the complex 
system of the shock wave and vortex sheet interactions which develop from the blade 
leading edges during the turning process. The flow system in a supersonic inlet 
channel must swallow the initial shocks generated by the leading edges of the blades. 
The bow shocks will eventually interact and a vortex sheet (or slipstream) or new 
shocks and/or expansion waves will be generated which will cause various combinations 
of shock wave interactions further downstream (Figure 2). To determine the true 
performances of a supersonic cascade, these interactions must be accurately calcu- 
lated. 


Near the two surfaces of the channel the flow is also complicated by the exist- 
ence of the boundary layer. On the concave or compression side of the walls, the 
strong positive pressure gradient can cause boundary layer separation which greatly 
decreases the channel efficiency. Shock wave-boundary layer interaction occur when 
shocks impinge on the blade surfaces. The sudden pressure rise across the shocks 
can easily cause boundary layer separation or even boundary layer transition, Also, 
the fact that the boundary layer displaces the inviscid flow field gives rise to 
viscous-inviscid interaction effects. 


Another problem arises from the fact that the velocity ia discontinuous 
across a vortex sheet. In actual flows, this velocity discontinuity is modified 
by viscous mixing. The mixing problem, however, is beyond the scope of the 
present report. 


This report describes a FORTRAN language computer program, developed for the 
IBM 7090, for calculating inviscid and viscous flow fields within two-dimensional 
supersonic curved channels. The channel surface leading edges can be either sharp 
or blunted. The complete channel flow field is divided into three regions; each of 
which is calculated by an independently developed numerical procedure. These pro- 
cedures are chained together in the overall program to provide automatic calcula- 
tion. These three regions are the blunt body, supersonic, end boundary layer flow 
fields. The blunt body flow field, located ahead of a blunt leading edge is com- 
puted by the second order Pelotserkovskii method (Reference 1), the supersonic 
flow field is computed by the method of characteristics, (References 4, 5 and 6), 
and the boundary layer or the viscous flow field is computed by a finite difference 
procedure (References 2 & 3). Singularities that can be treated by the program 
are: interactions of shocks of opposite family; vortex sheet; interactions of 
shocks with a vortex sheet; the wall of the boundary layer; and viscous-inviscid 
interaction effects which are directly related to the boundary displacement. It 
is assumed that the inviscid flow field is supersonic throughout, except in the 
blunt leading edge region. 


Three sample calculations are presented. Mach 4 and Mach 6 cases were computed 
for a sharp leading-edge channel, and a Mach 4 case was computed for a blunt leading- 
edge channel. 


Il CHANNEL AND FLOW MODELS 


The flow channels as formed by a cascade are shown in Figure 1. The blades or 
the turning vanes being treated can have either sharp or blunt leading edges. 
Typical design values initially given for studies are the incident angle a, turning 
angle $8, and the incident Mach number My. In the numerical calculation, the 
channel and flow characteristics are defined on a cartesian coordinate system x and 
y, where x is parallel and y is perpendicular to the incident flow direction. 
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Figure 2 shows a schematic diagram of the channel flow model considered in 
the present computer program in which it is assumed that the inviscid flow field 
remains supersonic downstream of the leading edges. Two initial bow shocks will 
be generated by the two leading edges. For the cases with blunt leading edges, a 
blunt body procedure is required to calculate the initial portions of the bow 
shocks along with a special procedure to calculate the two starting supersonic 
characteristics. For sharp leading edges, the usual wedge flow solutions apply. 
The bow shocks will eventually intersect as they extend downstream. At the 
shock intersection point, a vortex sheet or slipstream will be created. When 
the shocks impinge on the blade surfaces, they will be deflected again into the 
flow field and interact with the vortex sheet. At the point of shock interaction 
with a vortex sheet, a new expansion (Prandtl-Meyer fan) or a compression (shock) 
wave will be created depending on the flow conditions. 


When the inviscid fiow solution has been determined, the boundary layers 
developed on the blade surfaces can be calculated. Shock wave-boundary layer 
interactions are treated by spreading the pressure jump (or velocity discontinu- 
ity) across the shock along the surface. The spreading distance is a function 
of the local boundary layer thicknesses. Viscous-inviscid interaction effects 
are treated by increasing the channel width locally by the boundary layer dis- 
placement distribution. 


III EQUATIONS AND NUMERICAL TECHNIQUE 


Three basic numerical procedures are used in constructing the overall com- 
puter program. These procedures are: the Belotserkovskii method for calculat- 
ing flows around the blunt leading edges, the method of characteristics for 
calculating supersonic flow fields, and a finite difference procedure for cal- 
culating boundary layers. These methods and their equations are described: below, 
Since the basic theory of the method of characteristics has been extensively 
described in the literature (References 4, 5, and 6), the basic inviscid flow 
equations have not been included. The calculation procedures, however, are sum- 
marized in the appendices. Since the methods described here have been described 
previously in the literature, the symbols used in the equations are essentially 
those of the originals. Since no attempt has been made to generalize the nomen- 
clature, a list of symbols is given for each of the methods. 


A, BLUNT BODY S LONS 


The flow about the blunted nose is determined by using the direct method 
developed by Belotserkovskii (Reference 1). This method is based on a general 
technique proposed by Dorodnitsyn (Reference 7) for solving a system of first 
order partial differential equations of mixed type which have two independent 
variables. The technique consists of reducing the partial differential equations 
to ordinary differential equations in one independent variable by Specifying the 
variation of the unknowns with the other independent variable. 


The ordinary differential equations which result are integrated numerically. 
Most of the boundary conditions are given at the beginning of the range of integra- 
tion but these must be supplemented by conditions that the solutions of one or 
more of the equations shall pass smoothly from the elliptic (subsonic) to the 
hyperbolic (supersonic) region. The conditions of smooth transition, together 
with the other boundary conditions, determine solutions of all unknowns uniquely 
(Reference 8). 
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Since the nose of the body is a circular cylinder, polar coordinates with 
the origin taken at the body geometric center are used. Figure 3 illustrates the 
geometry of the problem. Dimensionless variables are used whereby length is re- 
ferred to cylinder radius, velocities are referred to the limiting speed, and 
pressure and density are referred to corresponding stagnation values in the free 
stream. 


The fluid motion is governed by Bernoulli's equation, the equation of contin- 
uity, an equation obtained by combining the radial momentum equation with continuity, 
the condition of conservation of entropy along streamlines and the equation of state, 
The equations may be written: 


Continuity equation 


F ot _ (3a.1) 
yz (th) + 95 =0 


Modified momentum in the r-direction 
é os 3 
ee ———=> ap = a.2 
In equations 3a.1 and 3a.2 


H kp + pu- 


8 Ppuv 


g kp+pv 


h 


Tu 


t 


TV 
Bernoulli's equation 


> = ple w) anit 


Conservation of entropy 
¢ = pl” = (yp) (3a.4) 


where W is the stream function defined below, and # is a function of entropy. 
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FIGURE 3. GEOMETRY OF BODY AND SHOCK 
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The equation of continuity leads to the stream function defined by: 


oY _ 
oy = Feu 
(3a,5) 
ay _ 
er ~ Vv 
From Eqs (3a.3) and (3a.4): 
as (32.6) 
- a. 
p =r%y 
1 
a= © (3a.7) 
p = T¢# 
1 
where: T=(l- w’) yl 
From Eq (3a.5), on any line: 
ad e pw 55 - ru). (3a.8) 
The speed of sound is determined from Bernoulli's equation: 
sy Pet (l-w). (3a.9) 
The boundary conditions on the surface of the cylinder (r = ro = i) are: 
u=wW=0 
and: 
2 2 
4y [y-l 1 Ww (y-1) (3a.10) 
y a ——f-y ge, Gein! 2 
v2.1 =) we? lew 4y |. 
The equation of the shock wave is: 
ro o=7; (0) = 1+ €(6) (3a,11) 


where €(@)is the distance from the body surface to the shock wave along a line 
6 = constant. 
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The conservation conditions across the shock wave lead to the following relations. 
The velocity components behind the shock wave parallel and normal to the free stream 
velocity vector are: 


2 2 1 
Ww, = W, | 1-= sino - (3a,12) 
= e* ( rh x) 
1 1 
: ¥ me sinto) 71 


The radial and transverse velocity components are given by: 


u, = Wy Sin 6 ~ Ww, cos 6 


v; w, sin 6 + w, cos 6. (3a.13) 


The pressure and density are given by: 


y 2 
_ 4y Ayr | 8 2 (y21)* 
oa Oo Ee sgt] on 


1. wo 
1 2 
a rel ( : w2) y-I een eee (3a.15) 


1+(1- we ) cot“o°. 
Since ¢ = p/p” » one obtains: 


sooth Gal bs] af oan 


where: 


lew, 


The value of the stream function ¥ is determined from the condition of con- 
tinuity across the shock wave. Thus: 


1 
vey = Ww (1-w2)77F sedan, (32.17) 


From the geometric relationship on the shock wave, & =tano , one obtains 
the equation for €(6). 


ay = = (1 +€) cot (0+6) (3a.18) 


The approximating system of equations is obtained as follows. Introduce the 
new independent variable: 


t = — (3a.19) 


After changing independent variables from (r, 0) to(é,8) and integrating Eqs (3a.1) 
and (3a.2) along a line @ = constant from ¢=Q to é= Sy one gets: 


FI 
é = 

rh -r; tts / etdgé =0 (3a.20) 

° 

st 4 

3 = 

ee coat [ eg dé. (3a.21) 
O 0 


To obtain the approximating system of ordinary differential equations, the layer 
between the body and shock wave is divided into N equally spaced strips and the 
integrals in Eqs (3a.20) and (3a.21) are evaluated on successive boundaries of these 
strips. The unknown functions in the integrals are approximated by interpolation 


polynomials of degree Nin € taking the interpolation points to be the strip bound- 
aries: 
N 


{(¢,0) = 2, a, (6) 8%, (Sa.22) 


where the coefficients a,(6) are linear functions of the values of f on the strip 
boundaries. 


Writing Eqs (3a.20) and (3a.21) along each of the N-1 intermediate lines, 
é a Noi+t (y and 

i N ¥ on the shock wave and on the body are determined from 
the boundary conditions), and taking Eq (3a.18) into account, one obtains an 
approximating system of 4N-1 equations for the unknowns: 


€,y oO, Vor Uy Vys Vy» v, (i= 2, 8, ooo ND. 


In Reference 1 it was shown that the convergence of this method with increasing 
N is extremely rapid, so that N= 2 is sufficient for most purposes. This is termed 
the "second approximation" and is the one that is used here. 


In the second approximation, the system of equations is: 


ae =z =~ (1 + €) cot (¢+ 6) 


is 7 F, (3a,23) 
dv E 
0 0 
z (8a.24) 
qo “Do 
where: 9 
2c 
te) ' _ y-l 2 
EY y+1) 7, 6° Do = HT - Yo 
dv, 1 € 
ais Sy S,' = Uy , yi t,' a v1 av. a6. (3a.26) 
to % ~ yal : 
dv E 
2_ 2 (3a.27) 
qe “Do 
where: 
2 
Cc t,u 
2 2 ; 2 
E> on Tr) ( +—7- ss 
C9 
Yel 2 2 
Do=ss1 7 Yo + Far “2 
Yo = # (Wo) (3a.28) 


Quantities appearing in Eqs (3a.18), (3a.23-3a.28) are given by: 
2 ey SEP 1 


i 2 €' 4 aa 
S, = (3 S, - 485) - = | Hy + (1 +e) H (2+e) Hp | +8, Eo (3a,30) 
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where: 


1 


5, =i fe 8) +8 Hy 202+) Hse) Hy] + By +98, 


to’ = ty! + (4 to = ty = Sty) +2 [Ob sey my = + edt ] 


t.' 
= 1 e' 2+ 5 l+e 
to == +2£ (ty -t,) + ( € )o,-§ hy 


qyn . do do dé _ 
Yo= Vy 
df ( <1)” 
n¢ W = 
—r—> = 2 cot = 1 
(l+w) w - &¢=) 
ayy 2 pat 
“Ae = We (l- w,, a | €* sin 0 (1 +6) 08 6 | 
: u, V, sin 2 
D, -4 wz (1-w_2) 72 1 "1 
1 


1-w, 


2v 
1 
pf as Ga | yy OF a te m,)]| 
ae | 


V 
i 

G) #7; 2% ss aes | =) 4] 
1 


dw dw 
m, oer ead sin 6 --7> cos 6 


(3a.31) 


(3a,32) 


(3a.33) 


(3a,34) 


(3a.35) 


(3a.36) 


(3a.37) 


(3a.38) 


(3a.39) 


(3a.40) 


dw, dw ‘“ 
ny =--7> sin 6 ---~ cos 6 (3a.41) 


dw 2w 
a> =. =r sin 20 (3a.42) 
dw 2w (32.43) 
= cos 20 + ) a. 
a ( Mm sin’o 


The integration of the above system is carried out numerically and begins at the 
axis of symmetry 9=(Q, where g => v = Vo = Vo = Q while the values € (0) and uy(0) 
are unknown. 


Examination of the equations for vy_ and y, reveals that singular points exist 
where: 0 2 


2 

Zu 
2 yel 2 yal 2 
‘Spr Oe ee ee 


that is,where the transverse velocity components on the body and on the mid-way line 
become sonic. In order for a continuous solution to exist, one requires that: 


Eo =Q when D, =0 


(3a.44) 
Eo =0 when De = 0, 


Otherwise %S and Vo become infinite at the singular points and the solution is 
physically meaningless. Thus the two conditions of smoothness determine the proper 
values of €(0)and U,(0). Since one is dealing with a two point boundary value problen, 
iteration mist be Used €6 determine the correct values of €(0) and u,(0). Once this is 
done the numerical integration is contimed a short distance into the supersonic region 
where an initial value line is generated to start the characteristic procedure. 


The initial value line is a set of 20 equally spaced points along a line @ = 
constant between body and shock wave. Since only the flow quantities on the body, 
mid-way line and shock are known from the solution, additional quantities in between 
are obtained by fitting quadratics to u, y and @. All other flow properties may be 
determined from these three, 


B, METHOD OF CHARACTERISTICS 


The supersonic inviscid flow field is solved by the method of characteristics, 
with the calculation proceeding along left-running characteristics, Special treat- 
ment is required for singularities in the form of shocks and contact surfaces. In 
the following discussion the characteristics and shocks will be labeled as follows: 


Right running characteristics of the second family 
Left running characteristic of the second family 
© left running shock of the first family 
@ Left running characteristic of the first family 
@ Right running characteristic of the first family 


3 Left running shock of the second family 


Starting Region For Pointed Blades 


Initial values for the method of characteristics on pointed blades are obtained 
from the solution for the flow on a two dimensional wedge. The shock angle at the 
leading edge of the blade is determined by incident flow and leading edge angle. 

The procedure then eaters a body point procedure and erects a right running char- 
acteristic from the shock to the body. Using this body point a new shock point 
is calculated which is used to calculate a new field point. The order in which 
the points are calculated is shown in the sketch below: 
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Left running characteristics are built up until they intersect the shock wave. 
In this way the inviscid flow field behind the bow shocks can be calculated up to the 
point where the shock waves interact. 


Starting Region For Blunt Blades 


The determination of the subsonic flow region of blunt blades has been 
described in the previous section. The Belotserkovskii solution can be used for 
a short distance into the supersonic flow region to provide an initial value 
line for the method of characteristics solution. The blunt body solution pro- 
vides values at three points: at the shock, midway between the shock and body, 
and at the body. Twenty points are then generated on a straight line running 
from the shock to the body using quadratic interpolation, 


STARTING 


INITIAL VALUE CHARACTERISTIC 


LINE 
—— 


SONIC LINE 


mii 


The initial value line is chosen such that the minimum Mach number of the three 
points from the Belotserkovskii solution is 1.1 The characteristics network is then 
built up until a complete left running characteristic is formed. After this the 
calculation proceeds as before until a shock interaction is reached, 


Shock Interactions 


a. Interactions Of Shocks Of Opposite Families 


When two shocks of opposite family interact, two reflected shocks and a slipstream 
are produced: 


PRIMARY REFLECTED SHOCK 


SLIPSTREAM 


REFLECTED SHOCK 
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The solution for the two reflected shocks must be iterated to determine a 
common flow deflection angle ¢ behind each shock, 


A good firet approximation for this turning angle is the difference of the 
two flow turning angles across the primary shocks, The angle ¢ is altered in 
the iteration process until the pressures match across the two reflected shocks. 
A similar procedure is used to calculate each point along the slipstream. 


Conditions to be satisfied across the slipstream are that the pressure and 
flow direction match across it. Discontinuities exist, however, in temperature, 
velocity and stagnation pressure. 


b. Shock Interactions With The Wall 


The problem of calculating the interaction of a shock with a surface resembles 
that of calculating the interaction of two shocks of opposite families in that the 
surface corresponds to the slipstream. In the case of the interaction with a surface, 
however, the stream direction after the reflection is determined by the wall. The 
reflected shock is such that it matches the flow previously calculated downstream 
of the impinging shock with the flow calculated along left running characteristics 
after the reflected shock, 


c. Shock Interaction With A Slipstream 


A shock interacting with a slipstream will produce a transmitted shock and a 
reflected wave, either expansion or compression: 


TRANSMITTED SHOCK 
—_——_ 


SLIPSTREAM 


IMPINGING SHOCK ‘eo REFLECTED WAVE 


As in the case of the interaction of shocks of opposite families, an initial slip- 
stream direction is chosen which determines the angle (and strength) of the trans- 
mitted shock. The first approximation to the slipstream direction is taken as the 
flow direction after the impinging shock. The pressure matching condition is then 
imposed, and the flow direction is altered until this condition is satisfied. The 
flow direction immediately after the interaction point determines whether the re- 
flected wave is an expansion or a compression. 


Slipstream Calculation 


The flow direction and the static pressure must be the same across a slipstream. 
To calculate points along a slipstream, the flow direction must be iterated by apply- 
ing the compatibility equations along the left running characteristics from both 
sides, 


SLIPSTREAM 


RIGHT RUNNING 
CHARACTERISTICS LEFT RUNNING 


CHARACTERISTICS 


C. BOUNDARY LAYER SOLUTIONS 


Laminar boundary-layer solutions are obtained by a numerical method described 
in References (2 & 3). In this method, an explicit finite-difference technique is 
utilized to solve the exact boundary-layer equations for steady compressible flows. 
It is applicable to any smooth surface with arbitrary boundary conditions (pressure 
and temperature distributions, etc.). 


Since the procedure used in this program was originally a general three-~ 
dimensional program, the equations described here are the original boundary layer 
equations in three dimensions. The equations in general curvilinear coordinates 


€, » , § are: 
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where € and 9 represent a two-dimensional orthogonal curvilinear dimensionless 
coordinate system which is defined appropriately for the given body surface, and 
f is orthogonal to € and 9 and corresponds to the vertical distance. hh 
and ha are the metrical coefficients in theé, ” and § directions respectively, 
For two dimensional flows, all properties, of course, are independent of the 


" coordinate. 
The total enthalpy E in Eq (lc) can be eliminated by making use of the relation: 
—~ 2 +2 
ee ee (2) 
= ee 
Substituting this expression into (lc) and combining it with Eq (la), the 
energy equation becomes 


cat 
“p 


a a ee tee ee peas = om _ —arZ 
u et Vv @ w oe 1 u op V_ op 7) au 
ee ea ff += + —- 
n. °§ 7 On FF 8F mT (5 af + z ) —=,|(2) 
hy hy hg cA hy § he cA hg 

(2) | aa +(A4 (1e) 

é ~»nH 5. O¢ 
4 c)Ahs hg 
The coordinate ¢§ is defined as: 
t -(1-3)”? (3) 
U 


where u and U are the dimensional internal and external longitudinal velocity com- 
ponents respectively, Therefore, {§ varies, 0 to 1, from the outer edge of the 
boundary layer to the surface of the body. 


A shear coeff‘cient g is introduced, defined as: 
-_f# @ WU i 
hg 


The equation of state: 
p =ARt (5) 
is introduced to eliminate the density h , and the definition of { removes explicit 
dependence ony . As a result, the system of Eqs (1) is reduced to a system of four 
equations in the four unknowns v, t, » and W. 
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Here, A: * * * A18S are the boundary quantities that must be specified for the 
boundary layer integration. They are completely determined by the external inviscid 
flow condicions and the body surface geometry, They are evaluated in terms of the 
inviscid flow velocity vector distribution and the metrical coefficients in the & and 

” coordinates. Thus they are independent of the {§ coordinates. 


The dimensional physical quantities in Eqs (1) are replaced by their dimensionless 
equivalents which are defined as: 
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Subscript "o'" denotes reference quantities. Upper case letters denote boundary 
layer external quantities. 


ie is an arbitrary reference length; R is the local Reynolds number defined as 
UAL SH» O A singularity factor g* is incorporated in the quantities associ- 


ated with §, where by definition &* is identical to € in a segment adjoining a 


Stagnation or sharp leading point or line and is any convenient positive constant or 
function of & and # elsewhere. 


Finally w is eliminated by substituting its expression from Eq (6a) into the 
remainder of Eqs (6) yielding: 


2 
B, oF + Ba oz 7 + (Bs + Bs) = By ¢ 377 + Bg (7a) 
ot ot ov 2 
By $f + By 2 + (By + By) $y = By o 57a + By + Bo 23, 
(7b) 
2 


a2, 

Bi oF t+Bos —s By ¢ 5+ By + Bio or 7 

t t 
at a 

+ Big 3p + Bra oe i F | + Byg St + Bigop 


B, * * * Byg are differential coefficients defined as functions of the boundary 
conditions (surface geometry and inviscid flow properties), the transport properties 
of the fluid, and the non-dimensional quantities y, t, and ¢ : 
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For two dimensional flows the crossflow velocity V, and crossflow gradients 
8( )/an are equal to sero. 


In the three dimensional procedure, each boundary layer segment is replaced by a 
lattice of points with constant spacings a, b, c and integer indices 1, j, k in the 
&, 9, § directions respectively. In Eqs (7a), (7b), and (7c), first order deriva- 
tives appear with respect to all three coordinates, Es 9, and { ; second order 
derivatives appear only with respect to the { coordinate. In the central difference 
approximation, the first derivative difference equations are written as: 
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in Eqs (8a), (8b), and (8c) respectively. In a similar manner, the second derivative 
difference equation is written as: 
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This form of difference equation however, is usually numerically unstable for the 
parabolic type of equation (References 3, 9, 10). It was found in Reference 3 that 
numerical stability can be obtained if the second derivative difference Eq (10a) is 


rewritten as: 
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The requirement for this approximation, therefore, is that: 
a<<c, (12a) 


Other conditions which must also be satisfied in this procedure are: 


(1) a, b, c are small in general (12b) 
(2) a<<b (12c) 
(3) ah,lvisbh,u, (12d) 


Equations (7a), (7b), and (7c) therefore, can be written in the difference form 
by substituting Eqs (8a), (8b), (8c), and (10b) for the respective differentials. 
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where ¢ has been replaced by ck, and: 
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The numerical integration procedure is followed by a "leap-frogging" technique 
as described in Reference 3, 


Consider the point P; 4 , in the above sketch. The forward derivatives i,j,k 
in the i direction at Py j,, are calculated by Eqs (13a), (13b) and (13c) using values 
from the five nearest neighboring points shown on the ith and (i-1)th planes. The 
value of the new point ( 441,4,k 29 the (i+1l)th plane is calculated from the 
relation: 
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without calculating the value at P, j,k" 
$J 


In this procedure, therefore, if the indices i, j, k start from 0,0,0, the for- 
ward derivatives are calculated from Eqs (13a), (13b), and (13c) when (i+j+k) is odd, 
and the values are calculated by Eq (14) when (i+j+k) is even, 


IV PROGRAM CONSTRUCTION 


This program contains ten chains for the various phases of computation. Figure 18 
is an overall diagram showing the linkage of the system, The functions of each chain 
are as follows: 


CHAIN I 


Chain I reads in the first data card which contains a blade type indicator; free 
stream Mach number, pressure, and temperature; and the specific heat ratio . These 
data are stored on tape 2 for use by the entire program. If a blunt leading edge 
case is indicated, Chain VII is immediately called to calculate the flow field ahead 
of the blunt noses, If a sharp leading edge case is indicated, the wedge program is 
used to determine three initial field points for each surface as shown in the follow- 
ing sketch. 


Upper Blade 


Lower Blade 


The calculation of these points are controlled by Subroutine WEDGE. The first 
point (1) is calculated by Subroutine SHKL, and the other two Points (2) and (3) are 
calculated by a special shock-body subroutine (Subroutine SHOCK). The procedure of 
calculation used in these routines are described in Appendix A. 


The blade geometry data are then read from cards, and the surface data are 
fitted with a series of cubics using Subroutines CURFIT and CUBIC. The coefficients 
of these fitted equations are stored on tape 2 for later use in the progran. 


CHAIN VIL: 


Chain VII calculates the flow about a cylindrical body of unit radius by a 2nd 
order Belotserkovskii procedure. This chain determines the data along an initial 
value line in a radial direction from the body center of curvature to the shock, 
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FIGURE 16 PROGRAM CONSTRUCTION DIAGRAM 
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CHAIN VIII 


This chain uses the initial value line obtained in Chain VII and calculates a 


complete left running characteristic from body to shock, as shown in the following 
sketch. 


Starting Characteristic 


Initial Value 


The calculation starts at the shock, Characteristic I is first calculated, then II, 
etc., until a complete characteristic from the body to the shock is obtained. This 
starting characteristic can be used for the upper and lower blade surfaces by proper 
data transformation. 


The blade geometry data are then read from cards. The starting characteristics 
just calculated are dimensionalized with the given body leading edge radius. The 
blade geometry data are then curve fitted, merging smoothly to the cylindrical blunt 
leading edge sections. The coefficients of the body description equations are stored 
on tape 2 for later use in the program. 


CHAIN Il 


This chatn calculates the flow field of the upper and lower blade surfaces using 
the starting characteristics calculated in Chain I or Chain VIII, and continuing until 
the two bow shocks cross (Figure 2). 


CHAIN III: 

Chain III calculates the interaction of shocks of opposite families. The loca- 
tion of the interaction point is calculated by a curve fit method for intersecting 
two curves (Subroutine MEET), 

The properties behind each shock are calculated by Subroutine SHK1, and the slip- 


stream deflection and the properties behind each of the reflected shocks are calcu- 
lated by an iteration procedure which is discussed in Appendix A, 


Next, the first extensions of the reflected shock waves and slipstream are cal- 
culated, The sketch below shows the points being calculated in this chain. 
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CHAIN IV: 


This chain controls the remaining characteristic calculation, which includes the 
slipstream, shock-body interactions, and shock-slipstream interactions. 


CHAIN V: 


This chain reads the characteristic data from tapes 3 and 4 and calculates the 
surface velocity vs surface distance data. These data are then written on tape 9 
for the boundary layer calculation. 


CHAIN X:; 


Chain X is the boundary layer input program which calculates the "A" coefficients 
from the surface velocity data. These coefficients are the input data for the 
boundary layer integration program. 


CHAIN XI: 
Chain XI is the boundary layer integration program which calculates properties 


across the boundary layer (velocity, temperature and shear profiles). From these 
properties, the boundary layer thickness and skin friction, etc., are calculated. 


V__SAMPLE RESULTS AND DISCUSSION 


Three sample cases have been calculated. Phe flow in a sharp leading-edge 
channel was calculated for incident flows of Mach 4 and Mach 6 while only a 
Mach 4 case was calculated for a blunt leading-edge channel. The calculated 
flow geometry for these three cases are shown in Figures 4, 5, and 6. For the 
two sharp leading-edge cases, inviscid flow fields were calculated up to and 
including the shock-slipstream interaction point. For the blunt leading-edge 
case, however, the calculation could not be carried as far, due to the fact that 
adjacent characteristics overlapped, which indicates that envelope shocks are 
being formed in this region. As yet an envelope shock routine has not been 
developed, The reason for the envelope shock formation is due, of course, to 
the particular channel geometry selected. These results show that an envelope 
shock subroutine is needed to complete the present type of flow field calculation. 


The sharp leading-edge geometry selected for calculation has 20° and 15° 
leading-edge angles on the lower and upper surfaces, respectively. The initial 
shock strength, therefore, is higher on the lower surface. Also, because the 
flow along the lower surface is being compressed while the flow along the upper 
surface is being expanded, the lower shock strength will be even stronger at the 
shock interaction point. When two shocks of different strengths interact, the 
resulting flow is deflected toward the side generating the weaker shock. For 
the Mach 4 and Mach 6 cases, the flow deflection angle was found to be about 7° 
and 8°, respectively. The general flow pattern for the two cases are similar 
except that at the point of shock-slipstream interaction, a compression wave was 
created for the Mach 4 case, while for the Mach 6 case, an expansion wave was 
created. These waves were found to be extremely weak and can be neglected in 
the downstream flow field calculation. Although the calculations have not been 
carried fartherdownstream, the shock shapes indicated that interaction of shocks 
of the same family is likely to occur only slightly downstream of the last con- 
putations. This is due to the fact that the deflected shocks from the upper 
surface extend much farther downstream than those of the lower surface. 


Because the flows are being expanded near the upper surface and the flows 
are being compressed near the lower surface, large pressure gradients exist 
across the channel during the turning process. Figures 7 to 10 show the pressure 
and Mach number distribution along the upper and lower surfaces for the two sharp 
leading edge cases, and Figure 11 shows the pressure variation across the channel 
at various longitudinal stations for the M=6 cas only. 


A blunt leading-edge radius of 1.0 was selected for the blunt leading-edge 
channel geometry. The incident angle « for this channel is about 13.75° 
(Figure 6). For the selected body shape, the flow has turned to an angle of 
approximately 30° before two adjacent characteristics overlapped. 


The laminar boundary layers were calculated for the two sharp leading-edge 
cases, The boundary layer displacement thickness distributions are plotted in 
Figures 12 and 13, (note that the vertical scale has been highly magnified). The 
calculations show that for both cases the boundary layer on the lower surfaces 
had separated at about the same longitudinal station - X = 10.0 for the Mach 4 
case and X = 10.65 for the Mach 6 case, These boundary layer separation points 
are also indicated in Figures 4, 5, %Z and 8. Separation occurs at the point 
where the shear, which is proportional to the velocity derivative du'/de' at 
the wall, becomes zero, Some velocity and shear profiles along the lower ‘urfaces 
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FIGURE 6 INTERNAL FLOW FIELD OF A BLUNT LEADING EDGE CHANNEL 


O°y = “W ‘TANNVHD 3903-DNIGV21 d¥VHS V 40 SdVsuNS ARL ONOTV NOILNGIULSIG auNnssaud 


X - d0usIs{q 1 euTpPNITBu07 


< aunstid 


OFISY S1NSSOIg 130] WEOIIS 9824 02 2738IS 18907 


Layer 


; 
3 


OF IVY SANSse1g [BIO] BVIIIS 9023 03 2728IS [#907 


on Point 


bet 
a 
bet 
a 
@ 
~” 


Longitudinal Distance - X 


PRESSURE DISTRIBUTION ALONG THE SURFACES OF A SHARP LEADING-EDGE CHANNEL M,.. = 6.0 


FIGURE 8 


- ys ae alate Lert Art ERA AT SIME aes aimee = Ps 


0 


14 


aan 
| 
SN 
ae 
|| 


& 
6 8 
Longitudinal Distance ~- X 


a 
n” “ “ ~ = ond 


2 


© 


29QUNN YS edezZ2NS 


MACH NUMBER DISTRIBUTION ALONG THE SURFACES OF A SHARP LEADING-EDGE CHANNEL Ma = 4.0 


FIGURE 9 


«© 


16 


10 12 14 


8 
Longitudinal Distance - X 


6 


AOE 
tN Ey 


J9qUNY YOBy 99ezaNS 
% 


4 


2 


MACH NUMBER DISTRIBUTION ALONG THE SURFACES OF A SHARP LEADING-EDGE CHANNEL M,, = 6.0 


FIGURE 10 


0°9 = “~H ‘IZNNVHD 250g-SNIGVZT dUVHS V 80d 
SNOLLVIS ‘TVNIGALIONOT SNOLEVA IV NOLINGINLSIG ZANSSABZ ASMAASNVEL 


K- eoueisyg eszeasuezy 
cl A | IT Ol 6 g L 9 $ 


Ye 
oD, 
YM 


7 


7S 
io 


rs 


wee2r38dtI1S 


Tt @anold 


OFISY SANSS0Ig 1820] wEOIIG 982g OL 2F98IgE {8907 


TEE 
“COCA 
Sceaeen 
AL 

SaKe ) 
“COA 
“CCAD 
“COCA 


it ttt AL 
SRR 


= S $ 4 g ° 


® e e e ra . 


9 


? 


6 


Surface Distance - X’' 


4 


3 


BOUNDARY LAYER DISPLACEMENT THICKNESS DISTRIBUTION ON THE UPPER 
(EXPANSION) SURFACE OF A SHARP LEADING-EDGE CHANNEL 


1 


FIGURE 12 


a@ > S8OURITYL JUSMODETdSTQ 


Surface Distance - X 


BOUNDARY LAYER DISPLACEMENT THICKNESS DISTRIBUTION ON THE 
LOWER (COMPRESSION) SURFACE OF A SHARP LEADING-EDGE CHANNEL 


FIGURE 13 


» @ + SSeUX TUL QUewWSDETdSIq 


are plotted in Figures 14 to 17. The separation points indicated on the pressure 
distribution curves in Figures 7 and 8 show that the boundary layer separates even 
with a very weak adverse pressure gradient. This indicates that in order to main- 
tain a non-separated boundary layer, the pressure gradient mist be kept small. 

To do this, however, a very long channel is required in order to obtain a reason- 
ably high flow turning angle. A more practical solution which should be investi- 
gated in the future, is to use boundary layer suction so that the wall shear 
values can be controlled externally. 


In view of the results obtained with this program, the following conclusions 
can be drawn: 


1, To calculate the complete internal inviscid flow field with turning, 
two additional IBM subroutines are needed; (a) envelope shock, and 
(b) interactions of shocks of the same family. 


2. Boundary laver on the compression surface can be easily separated by 
the adverse pressure gradient. Boundary layer suction may be needed to 
maintain an attached boundary layer. 
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SHOCK POINT CALCULATION 


APPENDIX A 


SHOCK POINT CALCULATION 


Two major shock subroutines are used for most of the shock point calculations. 

One is for an oblique shock calculation which is used to initiate a shock wave at a 
point of discontinuity in flow direction. The other is for a general field shock cal- 
culation which is used to extend a shock from a previously calculated shock segment. 
There are special cases, however, such as the shock-boundary interaction and shock- 
slipstream interaction, where modifications of these two basic methods must be made. 
The two major shock routines are described in detail in this appendix and the differ- 
ences for the special cases are outlined. 


OBLIQUE SHOCK ROUTLNE 


In the following, quantities ahead of and behind the shock will be denoted by a 
subscript plus or minus sign, respectively. 


The oblique shock routine calculates a shock angle w from a given wedge or flow 
deflection angle 6 . With a known incident Mach number Miand the fluid specific 
heat ratio Y , w and § are related through 


2 
(y+ 1) M 
cot 6 = tan o Bs -1] A-1 


2 (M_*sin“w- 1) 


Eq. A-l is a transcedental equation in w , and an iterative solution mst be 
used. Using the Newton iteration procedure, # is calculated by 


6-6 
= i = 


w 


eee 


Where subscript i denotes the ith cycle of iteration, and 


2 
a5 _ 2cos" 3 M+ cos 2w+_1 


- M, (Y¥+ cos2w)+2 sin w 
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(M, * sin 2w- 2 cot w) M, *sin2u 
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M, (¥+ cos 2w) + 2 


Here wy44,) represents a new approximation to the shock angle. To start the 
calculation, the first approximation to the shock angle is taken to be the incident 
Mach angle, i.e., 


1/2 
Ww = tan”! : A-4 


° Mo 


The general oblique shock calculation procedure can be summarized as follows: 


1. Assume a shock angle Wo (Eq. A-4) 

2. Calculate a 68, by Eq. A-l 

3. Check the accuracy by comparing 5, with 8. [s-8:| > 1079 
is generally considered as satisfactory. 

4, When the accuracy check is satisfied, the solution is converged; otherwise, 
a new shock angle 44) is computed (Eq. A-2); 

5. Go ce Step (2) and repeat the calculation until convergence occurs in 
Step (4). 


Equation A-1 contains two roots in # ; one corresponds to a weak shock solution, 
and the other corresponds to a strong shock solution. Only the weak shock solution 
has meaning here. For a given set of M, and ¥ , the desired root of Eq. A-l will 
be in the range 


Wy 4u< Wmax . 


Where “may is the shock wave angle for which M_ = 1.0. Similarly, Smex is the 
deflection angle corresponding to Wys,y. 
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_1|M,~ sin(2w_._) - 2cot (w 
P a 1 + max m 


mas M, [7+ cos (20, ,)| +2 
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where 


me. +1 23. 
max = |! + )M_ 3-Y¥+ . 
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V(7+ 1) [rs DM, 4-2(8-7)M, 2 (rea) A-5b 


It is necessary therefore,to first calculate 6 max to make certain that 
& < Smax: Also, since d§/dw approaches infinity at Smay , 4 54/duy, 
in Ey. A-2 is arbitrarily restricted to a maximum value of 1,0 during iteration. 


With w calculated, flow properties immediately behind the shock can be deter- 
mined directly. These properties are: 


l. Static pressure 


2¥M_ sin” w- (Y-1) A-6 
eo eee 9 Saas 
2. Total pressure ‘ ; z 


(Y¥+ 1)M, +sin @ 
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GENERAL FIELD SHOCK CALCULATION 


The general field shock routine calculates extensions of an existing shock wave. 
The flow field immediately ahead of the shock can be either uniform or non-uniforn. 
The shock wave calculated in a uniform field is a bow wave and the shock wave caicu- 
lated in a non-uniform field is a secondary wave. 


Right Running Characteristics 
Ahead of the Shock 


Left Running 
Characteristics 


Right Running Characteristic 


Shock Wave Behind the Shock 


In the sketch above, a shock wave segment is to be inserted from 1 to 3. A 
double iteration procedure is required to calculate the shock angle w and the flow 
direction @ at 3. The general iteration procedure is as follows: 


First assume a linear extension of the shock wave at Point 1 to approximate 
the location of Point 3, which lies between 3A and 3B. Then by interpolation, 
using the data at 3A and 3B, and the oblique shock relation, the incident and 
the transmitted shock properties at 3 can be determined. The results are then 
checked with the compatibility equation along a right running characteristic 


from 3 to 4, Point 4 is also an interpolated point, which in this case lies 
between Points 1 and 2. 


The calculation for both cases (bow or secondary wave) is similar except that 
for the case of the bow wave, the locations of Points 3A and 3B have not been pre- 
viously established, but the flow properties are known (free stream conditions). 
So, for the case of the bow wave, Points 3A and 3B are determined by 


Xa4 = X, - £ cos (- H gy) 


¥,, = Yg + F sin(-H5,) 
B A-9 
Xeon = X, + > cos (- Ug) : 
¥.,=Y¥,+ & sin(-we,) 
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w3 and (8 +4)3 are related by the equation 

Ya°Y 

~— 1 - tan(w: 6 Tz 
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yz can be eliminated by writing 

Ya-Y 
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so that 


san Ya hy ane, Fe) ry = yg Ones he) a 
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X= A-10 


Eq. A-10 related X3 and (0,)3 for an assumed w3. Since all incident shock 
properties at 3 are interpolated, it is convenient to calculate the ratio RR1 = 
(X3 - X9,)/ (X3p - X3q) directly by the Newton iteration procedure. Thus, 


Yo, ~ Y X,,°-% Xen - X 
3A 1 3A 1 3B 3A 
{(RR1) = RR1 + ly. - RR1 +y——_ zr | YT. 
3B” “3A 3B ~ “3A/ “SB” ‘3A aa 
z(wo+w,) +6, +86 9p _8 
an [ 3° "1 lu SA. pri SA 
With RR1 calculated, the incident properties can be determined by 
Xy = Xen + RR1 (Xq,, - X, 4) 
Yg = Yoga + RR1 (Yo, > Yga) 
(P,)g = Pag + RRIU(Pop - Poy) 
(6,)3 = O34 + RRI(Osp - O54) 
A-13 


(H)g = Hg, + RRI(H 55 - gq) 
(R,)q = Rg, + RR1(Rgp - gq) 


Also, the transmitted shock properties can be determined from the following 
equations: 
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Before a new value for “3 is calculated, however, the currently calculated 
pressure (P-)3 is compared with the previous value to determine if the solution had 
converged, For this 


|e, a) ~ Ps 41) | <-002 Ps a) 


is considered satisfactory. To calculate Point 4, the Newton iteration procedure 
is used. Thus 


- x, \ [y, - ¥ - 
rama) «naa - (q2¥"] pox + (naa 3) 
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where 


With RR2 calculated, properties at 4 are determined as 


4 > X, + RR2 (xX, - X,) 


+ RR2 (Yo - Y,) 
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1 * RR2 (Py - P,) 


so that a shock angle correction Aw can be calculated by 


Peg So), * awagy wag (35) 


where 


(M,)q" sin” we ; sin 2 a, 


(5!) = (y+ 1) sin? 8, 


[35 sin’ Wg - ] sin 2 W, 


The new approximation to the shock angle, therefore, is determined by the 
equation 


3 (141) = gq) t4o A-18 


The process now is to recalculate the position of Point 3 using a new average 
value for flow properties. 


This iteration procedure continues until the pressure values in two successive 
iterations satisfy the criterion given above. 


The procedure for the field shock calculation is summarized as follows: 


1, Assume a shock angle w3. The initial approximation is taken as Ww). 


5 


2. Calculate RRL by iterating Eq. Al2 until RR1 - RR1 4.1) <10 


3. Calculate the incident and transmitted shock properties at 3 and determine 
if 
|(P_)s a) 7 (Pg a+1)| <-001 (P_)3 (i), 


4, If the above condition is satisfied the solution is converged; otherwise, 
calculate RR2 by iterating Eq. Al5 until 


-5 


5. Calculate a new w ;,), by Eq.'s Al8 and Al7, 
6. Go back to (1) and repeat the calculation, 


SPECIAL SHOCK-BODY, SHOCK-SLIPSTREAM CALCULATION 


In the previous section, the procedure for calculating a field shock was de- 
scribed, in which certain specific data are required for the computation. There are 
cases, however, that the same data are not available; these cases being the first 
extension of a shock wave which emanated from the body and the first extension of a 
shock which has arisen due to either the intersection of two shocks or the inter- 
action of a shock and a slipstream (see sketches below). 
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In these calculations, the data at Point 2 shown in the sketch of the previous 
section are not available, A body or slipstream Point 5 must be first determined 
by calculating a right running characteristic from Points 3 to 5 so that the values 
at Point 4 can be interpolated between Points 1 and 5. The body point and slip- 
stream point calculation procedures are described in Appendix C, 


INTERSECTION OF SHOCK OF THE OPPOSITE FAMILY 


The location of the shock interaction of opposite family is calculated by a 
special curve fit procedure. This procedure utilizes the coordinate points of the 
two shocks to determine the location and the respective slopes at the point of shock 
interaction. From these slopes the shock wave angles and the incident flow properties, 
the transmitted shock properties (M-, P-, § , and R- at Points 1 and 2), can be 
directly calculated by the oblique shock relationships given in section A-l. 


SHOCK-BODY INTERACTION Characteristic 
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This routine calculates a shock point at the body surface P3 from Points Py 
and P) in the above sketch. Here, the flow properties in regions (1) and (2) are 
calculated in a manner similar to the field shock procedure and the flow properties 
at region (3) are calculated by the oblique shock procedure. The incident properties 
in region (1), however, must be determined by calculating a characteristic from Py, 
to Ps (right runner in the sketch), The shock strength between regions (2) and (3) 


must be adjusted so that the resulting flow in region (3) is parallel to the solid 
surface, 


SHOCK-SLIPSTREAM_ INTERACTION 


Left Running (Transmitted) Shock 
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Reflected Wave 


Left Running Shock 


This routine calculates the reflected wave, the transmitted shock, and the trans- 
mitted slipstream resulting from a shock-slipstream interaction. The location of the 
interaction point is first determined in a manner similar to that of the shock-shock 
interaction procedure described in Appendix A-4, The properties at Points 1 and 2 
can be determined by interpolating the values along the incident slipstream which had 
been calculated beyond the point of interaction, Now, the procedure is to first as- 
sume a transmitted slipstream angle, and from this angle, the properties at Point 4 
can be calculated by the oblique shock routine (Appendix A-l) and the properties at 
Point 3 can be calculated by the field shock routine (Appendix A-2). The reflected 
wave can be an expansion or compression type depending on the change in flow direction 
from Points 3 to 5. The flow angle at Point 5, of course, is the assumed transmitted 
slipstream angle. The properties at Point 5 are determined by the expansion routine 
(Appendix D) when the reflected wave is an expansion wave; when it is a compression 
wave, the oblique shock routine is used. The condition to be satisfied here is that 
the static pressure across the slipstream (Points 4 and 5) must be equal. It is 
seen that for each assumed transmitted slipstream angle two values of pressure (Py, 
and Ps,) are generally obtained. The procedure, therefore, is to iterate on the 
transmitted slipstream angle until the pressure at Points 4 and 5 are equal. 
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APPENDIX B 
FIELD POINT CALCULATION 


Field points are points at which all flow properties are continuous. The method 
of calculation is to determine the point of intersection of two characteristics of 
opposite families passing through two known base points, and compute the flow pro- 
perties at this new point through the use of the compatibility equations, 


Streamline 


cil 
Left Running Characteristics 


Right Running Characteristics 


Streamline =e 3 


The above sketches show the two base Points (1 and 2) and the point to be 
determined (3). Point 4 is the intersection of the upstream portion of the strean- 
line passing through Point 3 and a straight line connecting Points 1 and 2. Since 
the total pressure is constant along a streamline, the total pressure at Point 3 
can be determined by interpolating between 1 and 2 to obtain the value at 4. The 


location of Point 3 is found by solving the equations of the characteristics passing 
through Points 1 and 2; 


®t 


1? 


i 
Ne 
tan(dsu)g-5 0s tan( tu y 


3 = 
tan(6twT F tan(6+u)5 3 


Yo - Y 
X3- X1+ EAT B-2 


where the average values (2,3,1,3) are first assumed to be the values at Points 1 
and 2. Now, t.e compatibility equations may be solved in their finite difference 


form for P3 and 93. 
Po + Ps 03 - 62 
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In order to compute #3, the stagnation pressure ratio (R) at Point 3, is 
needed, and is computed by first locating Point 4. Assuming R varies linearly 
between Points 1 and 2, R3z is determined by 


me ie Whar 
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The only unknown in this equation is X, which is calculated by 


X, =X, + AX 


ppail (X, - X})(tan 6 + tan 6,) aero 2(Y, - Y,) 
a) tan 6, - tan 6, X 3 tan 6, - tan 6, 


) (X, = Xs) (tan 94 - tan 90) - 2 (Yo = Y,)] 


These equations are derived from the following relations: 


tan 6, - tan ds KX, - X 


tan@,-tan@, X,-X, 


Y,-Y 


3 1 


which assume a linear variation of the flow Properties between Points 1 and 2. The 
root being used in equation B-3 is the one which makes | E a minimun, 
y~ = 
1 


Finally, the Mach angle at Point 3 ( My ) is calculated from 


1-7] 5 
y 


=] 1 


2 
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This completes one cycle of the calculation of_Point 3. These new flow pro- 
perties are used to determine new average values (1,3, 2,3) which are substituted 
back into the geometric equations (B-1 and B-2) to recalculate a new location of 
Point 3. The calculation is repeated until 


Paci) ~ P3i-t) 3(i) 


at which time the calculation of the field point is complete. 
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BOUNDARY POINT CALCULATION 


Boundary points considered here are those being calculated along the solid sur- 
faces and slipstreams. 


SOLID BOUNDARY 


Left Running Characteristic 


Right Running Characteristic 


Body 
~ 3 


In the above sketch, a body Point 3 is to be determined by calculating along 
a characteristic from 2. Since the total pressure is constant along a solid surface, 
the total pressure at 3 is obtained from the known value at l. 


The procedure is to first approximate the location of Point 3. The body curvature 
is defined locally by a cubic in the form 


Y = AX? 4 BX? +CX+D, C-1 


and the characteristic equation between 2 and 3 is given as 


Y.- 


3 2 C-2 
Xe = X 9 
Combining, we have 
3 2 
f(X,) = AX, + BX, + CX, + D- Yg - (X,-X,) tan( 6+ uss C-3 


from which X3 can be determined by a Newton iteration procedure. With X3 computed, 
the properties (@, P, #& ) at 3 can be determined directly by 


g,=tan”’ (Y4) = tan” (3AXg + 2BX, + C) 
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With these properties determined, a new set of average properties (2, 3) can now 
be calculated, and hence, a new approximate location of X43 can be determined, The 
calculation is repeated until 
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In the sketch on the previous page, Points 3 and 3' on the slipstream are to 
be calculated from Points 2 and 2' respectively. The procedure for calculating points 
on slipstream is similar to that for a solid surface except that here, the boundary 
location is not given and must be determined by matching the condition calculated 
from both sides of the slipstream. The conditions to be satisfied along a stream- 
line are that the flow direction and static pressure must be equal. 


Each segment of the streamline is assumed to be a straight line. Thus, by assum- 
ing a slipstream direction (from Point 1 to Point 3') the properties at Points 3 and 
3' can be calculated by a procedure similar to that of a solid boundary. In order to 
compare the pressure at the same location along the slipstream, however, a P3" is 
interpolated between Points 1' and 3' which should be equal to P3. In this way, by 
assuming a series of slipstream directions, a table of P3 and P3"' vs the slipstream 
angle can be established. The solution of the slipstream angle will be the value at 
which P3 = P3", This value can be calculated by the special curve fit procedure 
(Subroutine MEET) which determines the intersection point of two given curves. 
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APPENDIX D 
EXPANSION WAVE CALCULATION 


The configuration considered in the present program is shown in the sketch 


below where a discontinuity in flow direction is caused by a shock slipstream inter- 
action. 


-Expansion 


The calculation here is to determine the flow properties at 2 for a given 
flow deflection angle § which is equal to the increase in the Prandtl-Meyer angle 
v(M) from 1 to 2. The Prandtl-Meyer angle is defined as 


y(M) = y aa | tan = (M2-1) - tan! y M1 D-1 


From a known Mach number at 1 (M)), ¥) (M,) can be directly calculated. To 
determine M> for a given deflection angle, however, an iteration mst be carried 
out. The flow deflection §, is related to the Prandtl-Meyer angle by 


$= v (M,) = vi (M,) D-2 


The method of calculation is to first establish a table of M vs vy (Eq. D-1) 
and by interpolating this table, Mp can be determined directly froma V, calcu- 
lated by Eq. D-2. Thus, from Mo, the properties at 2 can be calculated from 


= y-1 
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Ho = tan i 
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where Ra = R,, since the flow expansion from 1 to 2 is isentropic. 
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AUXILIARY SUBROUTINES 


SUBROUTINES CURFIT ~ 
Subroutine CURFIT computes the coefficients of the curve-fit function 
y= Ke 


for a given table of y vs x such that the first and second derivations [fy'(x), £4 "(x)] 
are continuous within the range of the given table. The function £4(x) is defined as 
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which is a cubic expression defined by the points at i and i+l, and the two coef- 
ficients (a, ,b,) for each of the given if intervals (see sketch below). 
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The coefficients ay and by are determined by setting the first and second de- 
rivatives of the interpolation formulae for the ith and (i+l) th intervals equal at 
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Equations E-4 and E-5 apply for the i + 2, n-2 intervals, and hence, represent 
a system of 2n-4 equations in 2n-2 unknowns. Two additional equations, there- 
fore, are needed, and they can be obtained by specifying either the first or second 
derivatives at both ends of the table. These equations are 
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Equations E-4, E-5, E-6 and E-7 represent a tri-diagonal matrix system of 
the simultaneous linear equations in a, and b;. Hence, the solution (a;, by) can 
be obtained by a special matrix inversion procedure which is done in CURFIT, 


SUBROUTINE CURVE 


Subroutine CURVE utilizes the curve fit coefficient computed in subroutine 
CURFIT to calculate the quantities y(x) and y'(x) for a given x. The subroutine 
will also perform straight line extrapolations beyond both ends of the table. 
SUBROUTINE CUBIC 


Subroutine CUBIC recalculates coefficients computed by subroutine CURFIT to 
those corresponding to equations of the form 


y(x) = Ax® + Bx? + Cx + D 


A combination of straight lines and cubic curves can be also fitte*: smoothly 
within a given table. 


SUBROUTINE MEET 


Subroutine MEET determines the intersection point of two given curves which 
are specified as two independent tables y),(x) and y(x) 


y, (x) 


Intersection Point 


By using subroutine CURFIT, a table of x vs Ay(x) can be first established, 
where 


4 y(x) = yo(x) - y,(x) 


This table can then be curve fitted to obtain x, of the intersection point where 
Ay-o. With Xo determined, therefore, the , and the derivatives y') (x,.) and 
Y¥' 9 (x9) at the intersection point can be again determined by subroutine CURFIT. 
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APPENDIX F 


DESCRIPTION OF DATA INPUT 
AND OUTPUT 


This program requires the following input data: 


(1) 
(2) 
(3) 
(4) 
(5) 
(6) 
(7) 


A title or job description card, 

Blade type indicator (sharp or blunt leading edges). 

Mo - free stream Mach number. 

Pe - free stream pressure (atmospheres). 

Teo - free stream temperature (°K). 

Y  - specific heat ratio. 

Channel geometry - in numerical form, Two X vs Y tables must be provided, 
one for the upper (expansion) surface and the other for the lower (com- 
pression) surface. For blunt leading edge cases, the leading edge radii, 
the location of the radius centers, and the tangent point angles ( @y and 
Q,, in degrees) are also required. The formats and the ordering of the 
data are described in Table F-1, 


In preparing the input data, the following assumptions and rules should be 


observed: 


1) 
2) 
3) 
4) 


5) 


6) 


Ny, Ny cannot be greater than 20 for a sharp leading edge case and cannot 

be greater than 18 for a blunt leading edge case. 

Ry = Ry. 

61> by > 65°. 

To be compatible with the present boundary layer program, y should be set 
equal to 1.4. 

For a sharp leading edge case, straight lines are fitted for the first and 
the last segments in the two X vs Y tables (between Points 1 and 2 and 
between Points N-l and N). If the surfaces are curvilinear in those regions, 
these points should be closely tabulated so that the straight line approxima- 
tion is valid. 

For a blunt leading edge case, the first tabulated point on each surface is 
the tangent point, and the first interval should be greater than one leading 
edge radius. This is required to fit a smooth curve in that region. A key 
punch form and two examples of data input are included in the following: 
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OUTPUT 


I INVISCID PROGRAM 
(1) Job title 
(2) Moo 
Poo (atm) 
Too (°K) 


(3) Body Description 


x ¥ 


|> 
Tes 
or 
Is 


where (X, Y) - input points 
(A, B, C, D) - coefficients of cubics 


(4) Characteristic Output (along a characteristic) 


Xx 
P - P/PL OD 


THETA - Flow Direction (radians) 
ZMU Ss ~~ Mach Angle (radians) 


Indicator - defines point type - shock, body, vortex, or field 
(field points are denoted by ***) 
W - shock angle (radians) 


Y P THETA ZMU R LNDICATOR 
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Il BOUNDARY DATA FOR BOUNDARY LAYER CALCULATION 


(1) Boundary Layer Input Data - Upper surface 


T2T1l 
P2P1 
V2U1 
P3P2 


across the primary shock 


pressure ratio at leading edge 


velocity ratio 
pressure ratio across shock body interaction point, 


temperature io | 


Hop okt 


Xx ¥ Surface Distance (X') Velocity (U'/Q»9') 


(2) Boundary Layer Input Data - Lower surface 


T2TL = temperature ratio across primary shock 
P2P1 = pressure ratio at leading edge 


U2U1 = velocity ratio 
P3P2 = pressure ratio across shock body interaction point 
x Y Surface Distance (X') Velocity (U'/Q,') 


III BOUNDARY LAYER PROGRAM 


(1) 


Shock 
T(1) * T(2) 
M(1) P(1) * P(2) 
UCL) * U( 2) 
* 
* 
* 
a * a at eel 
free stream properties Properties across 
primary shock 
where 
T - temperature (°K) 
P - pressure (atm) 
U - velocity (ft/sec) 


(2) RE/FT - Reynolds number per ft. 


(3) LO) - Lo! 
SOLAR HEAT (set to zero) 
EMISSIVITY OF SURFACE (set to zero) 
DEL 1, DEL 2 (Surface distance over which pressure jump is being 
spread fore and aft of a shock boundary layer 
interaction point) 


(4) INVISCID DATA 


S(surface distance) * (U! ‘ 


Shock boundary layer interaction comment 


(5) 


(6) 


(7) 


FIX DATA 


Ss (é 4) Uk (U'/Q5') 


| ars rams | ee | 
where 
Lg - integration step counter in the longitudinal and 
transverse directions 
V - inviscid transverse velocity component 
T - non-dimensional temperature defined in List of Symbols 
A - coefficients (Boundary Layer Equations) 
X1-SQ( & *) — 
CP ( Cp) - pressure coefficient = a 


a as 


U* (U'/Qo) - velocity ratio 


tT 
F(O)* 


- inviscid temperature (°K) 
- suction velocity (set to zero) 


I, J - defined in (6) 
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vertical step counter 

transverse velocity component (v) 
non-dimensional temperature (tr) defined in List of Symbols 
non-dimensional shear parameters (¢) (Boundary Layer Equations) 
non-dimensional viscosity (#) 
longitudinal velocity component (parallel to Q) 
transverse velocity component (perpendicular to Q) 
vertical distance 


\* (displacement thickness) 
&** (momentum thickness) 
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